!> 蛙跳法计算量相对较少，二阶精度，与 wp2 相当，可用于力学求解
module seakeeping_leapfrog

    use seakeeping_kinds, only: wp => sk_real_kind
    implicit none

    private
    public :: leapfrog_type

    abstract interface
        !> 加速度函数
        subroutine func(t, x, v, a)
            import wp
            real(kind=wp), intent(in) :: t                  !! 时间
            real(kind=wp), intent(in), dimension(:) :: x    !! 位置
            real(kind=wp), intent(in), dimension(:) :: v    !! 速度
            real(kind=wp), intent(out), dimension(:) :: a   !! 加速度
        end subroutine func
    end interface

    !> 蛙跳积分法
    type leapfrog_type
        procedure(func), pointer, nopass, private :: fcn    !! 加速度函数
        real(kind=wp), allocatable, public :: v(:)          !! 真实速度
    contains
        procedure :: integrate, step, destroy
    end type leapfrog_type

contains

    !> 配置蛙跳法
    subroutine integrate(self, fcn, t0, x0, v0, a0, dt, n)
        class(leapfrog_type), intent(inout) :: self
        procedure(func) :: fcn                              !! 加速度函数
        real(kind=wp), intent(in) :: t0                     !! 初始时间
        integer, intent(in) :: n                            !! 长度
        real(kind=wp), intent(in), dimension(*) :: x0       !! 初始位置
        real(kind=wp), intent(inout), dimension(*) :: v0    !! 初始速度
        real(kind=wp), intent(out), dimension(*) :: a0      !! 初始加速度
        real(kind=wp), intent(in) :: dt                     !! 时间步长

        self%fcn => fcn
        allocate (self%v(n), source=v0(:n))

        call self%fcn(t0, x0(:n), v0(:n), a0(:n))
        v0(:n) = v0(:n) + a0(:n)*dt/2

    end subroutine integrate

    !> 蛙跳法，计算一定时间内的位置、速度、加速度
    subroutine step(self, t, x, v, a, dt, m, n)
        class(leapfrog_type), intent(inout) :: self
        integer, intent(in) :: n                                !! 长度
        real(kind=wp), intent(inout) :: t                       !! 时间
        real(kind=wp), intent(inout), dimension(*) :: x, v, a   !! 位置、速度、加速度
        real(kind=wp), intent(in) :: dt                         !! 时间步长
        integer, intent(in) :: m                                !! 积分次数
        integer :: i

        do i = 1, m
            t = t + dt  ! 长时间模拟采用累加，不适用索引计数
            x(:n) = x(:n) + v(:n)*dt
            self%v(:n) = v(:n) + a(:n)*dt/2
            call self%fcn(t, x(:n), self%v, a(:n))
            v(:n) = v(:n) + a(:n)*dt
        end do

    end subroutine step

    !> 销毁蛙跳法
    subroutine destroy(self)
        class(leapfrog_type), intent(inout) :: self

        deallocate (self%v)
        nullify (self%fcn)

    end subroutine destroy

end module seakeeping_leapfrog
